/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#ifndef __BINOMIAL_H__
#define __BINOMIAL_H__

#include <stdint.h>

/* http://en.wikipedia.org/wiki/Combinatorial_number_system */

static inline unsigned int 
binomial (unsigned int n,
          unsigned int m)
{
  if (m > n)
    {
      return 0;
    }

  uint64_t numerator = 1;
  uint64_t denominator = 1;

  for (uint64_t k = n; k + m > n; --k)
    {
      numerator *= k;
    }

  for (uint64_t k = 2; k < m; ++k)
    {
      denominator *= k;
    }

  return numerator / denominator;
}

static inline unsigned int
set_to_index (const unsigned int* labels,   /* in descending sorted order */
              unsigned int        n_labels, /* length of labels vector */
              unsigned int        l)        /* total number of labels */
{
  unsigned int offset = 0;

  /*
   * does this make sense?  according to Mathematica
   * 
   * 1 + Sum[Binomial[n - m - 1, k - m], { m, 0, k - 1 }] == Binomial[n, k] 
   * 
   * is true.  That means the next index after the largest possible value
   * input to the second loop is the offset increment of the first loop.
   *
   * If labels is not sorted descending then this is fxored.
   */

  for (unsigned int m = 0; m < n_labels; ++m)
    {
      offset += binomial (l, m);
    }

  for (unsigned int m = 0; m < n_labels; ++m)
    {
      offset += binomial (labels[m], n_labels - m);
    }

  return offset;
}

static inline unsigned int
index_to_set (unsigned int  index,
              unsigned int* labels,
              unsigned int  l)
{
  unsigned int offset = 0;
  unsigned int n_labels = 0;

  for (unsigned int next_offset = offset;
       next_offset < index;
       ++n_labels, offset = next_offset, next_offset += binomial (l, n_labels))
    {
    }

  for (unsigned int m = 0; m < n_labels; ++m)
    {
      unsigned int labels_m = 0;

      for (unsigned int next_offset = offset;
           next_offset < index;
           next_offset = offset + binomial (labels_m, n_labels - m))
        {
          ++labels_m;
        }

      labels[m] = labels_m - 1;
      offset += binomial (labels[m], n_labels - m);
    }

  return n_labels;
}

#endif /* __BINOMIAL_H__ */
